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ABSTRACT 

The supermassive black hole of M87 is one of the most massive black holes known and has been the sub- 
ject of several stellar and gas-dynamical mass measurements; however the most recent revision to the stellar- 
dynamical black hole mass measurement is a factor of about two larger than the previous gas-dynamical deter- 
minations. Here, we apply comprehensive gas-dynamical models that include the propagation of emission-line 
profiles through the telescope and spectrograph optics to new Space Telescope Imaging Spectrograph observa- 
tions from the Hubble Space Telescope. Unlike the previous gas-dynamical studies of M87, we map out the 
complete kinematic structure of the emission-line disk within ~40 pc from the nucleus, and find that a small 
amount of velocity dispersion internal to the gas disk is required to match the observed line widths. We examine 
a scenario in which the intrinsic velocity dispersion provides dynamical support to the disk, and determine that 
the inferred black hole mass increases by only 6%. Incorporating this effect into the error budget, we ultimately 
measure a mass of Mbh = (S.S;!;^^) x 10^ Mq (68% confidence). Our gas-dynamical black hole mass continues 
to differ from the most recent stellar-dynamical mass by a factor of two, underscoring the need for carrying out 
more cross-checks between the two main black hole mass measurement methods. 

Subject headings: galaxies: active - galaxies: individual (M87, NGC 4486) - galaxies: kinematics and dynam- 
ics - galaxies: nuclei 



1. INTRODUCTION 

Supermassive black holes are thought to play a fundamen- 
tal role in the growth and evolution of galaxies (; Silk & ReesI 
[l998; Di Matteo et al. 2005). This idea is supported by the 
empirical relationships between the mass of the black hole, 
whose gravitational influence is limited to a very small region 
at the center of the galaxy, and the large-scale properties of 
the host galaxy, such as the bulge stellar velocity dispersion 
(g^; e.g., iF errarese & Merritt 2000; Gebhardt et al. 2000bl 
iTremaineet al. 2002; Giiltekin et al. 2009) and the luminos- 
ity (Lbui; e. g.. Dressier 1989; Kormendv & Richstone 19951 
iMarconi & Hunt 2003; Sani et al. 2011). Other interpreta- 
tions have also been put forth to explain the black hole mass - 
bulge correlations. For example, the relations may simply be 
the result of random mergers without the need for additional 
physical processes hk e active galactic nuclei (AGN) feedback 
(iPeng II2007I; iJahnke & Macci6ll201 lb . Our understanding of 
the interplay between black holes and their host galaxies is 
highly incomplete. Thus, accurately characterizing the shape 
and cosmic scatter of the relations is crucial for distinguishing 
between the various theoretical interpretations, as well as for 
determining the black hole mass function and space density 
of supermassive black holes. 

The black hole mass - host galaxy correlations are com- 
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posed of 87 black hole mass (Mbh) measurements, which are 
most often mad e by modeling the mo tions of stars or nu- 
clear gas disks dKormendv & Hoi 120131) . In spite of the im- 
pressive number of observations, the local black hole mass 
census is still incomplete, notably for low and high-mass 
black holes, and there remain major open questions. At the 
high-mass end, the slope, intrinsic scatter, and possibly the 
functional form of the Mbh - host galaxy r elations are not 
well constrained (Gultekin et al. 2009; Graha m & Scottll20 131: 
iMcConnell & Maii2013i) . Properly quantifying the scatter at 
the high-mass end of the Mbh correlations may be especially 
helpful as non-causal interpretations of the scaling relations 
predict a decre ase in the scatter at hig her black hole and 
galaxy masses (iJahnke & Macciol 1201 li) . Furthermore, the 
current Mbh - cr* and Mbh - ^bui relationships make strongly 
divergent predictions for the black hole mas ses in the most lu- 
minous and the highest-dispersion galaxies (|Lauer et al.l2007t 
iBernardi et al. 2007), leading to questions about which of the 
two correlations is more fundamental. At least one of the 
scaling relations must be wrong at the high-mass end, but 
currently it is not known which or how to reconcile the dif- 
ference. In addition, there have been very few consistency 
checks between the stellar and gas-dynamical methods, there- 
fore it is unclear how much of the intrinsic scatter is due to 
inconsistencies between stellar and gas-dynamical measure- 
ments, or if there a systematic difference between the masses 
derived using the two methods. 

Only direct comparisons between the stellar and gas- 
dynamical methods within the same galaxy can address some 
of these questions, but presently such tests have only been 
attempted on a very small number of galaxies. Stellar- 
dynamical black hole masses have been made for both IC 
1459 and NGC 3379, but the ionized gas turned out to 
be highly disturbed and a useful gas-dynamical Mbh mea- 
surement could not be obtained (Verdoes Kleijn et al. 20(m 
lGebhardtetani2000at IC^peUari et al.ll2002bt " lShapiro et all 



l2006t Ivan den Bosch & de Zeeu wll20To |). In NGC 4 151, the 
gas-dynamical determination ([Hicks & Malk:anll2008l) is gen- 
erally consistent with the stellar-dynamical mass, but the au- 
thors label their stellar-dynamical result as tenta tive because 
they w ere unable to find single best fitting model ( Onken et alJ 
l2007h . Also , the gas-dynamical m easurement for the black 
hole in M8 1 JDevereux et al.l2003l) agrees with a p revious, but 
preliminary, stellar-dynamical determination dBower et alj 
lapOO). This leaves six galaxies (NGC 3227, NGC 3998, 
NGC 4258, NGC 4335, M87, Cen A) for which the two types 
of black hole mass measurements can be meaningfully com- 
pared. The stellar and gas Mbh measurem ents for NGC 3227 , 
NGC 4258, and Cen A are consist ent (iDav ies et al.l 120061: 
Ristorini et al. 2007; Neumaver et ani2007t iH icks & Malkan 
2OO81; iSiopis et al.1 i2009c iCappellari et all 12009). while the 



stellar-dynamical mass exceeds the gas-dynamical determi- 
nation b y a factor of ^2 - 5 for the r emaining three 
galaxies (iMacchettoetn] 119971; iVerdoesKleiin et al.i 120021: 



de Francesco et all 120061; iGebhardt et al.l 1201 It IWalsh et alJ 



In particular, the supermassive black hole in giant el- 
liptical galaxy M87 has been the subject of numerous 
m ass determinat i ons, b egin ning with the pioneering work 
of ISargent et all (119781) and lYoung et al.] (119781) . who sug- 
gested the presence of a ~(3 - 5) x 10^ M© black hole 
based on the surface brightness and stellar velocity disper- 
sion profiles. There were concerns associated with these and 
other early stellar-dy namical black hole mas s measurements 
(see review by Kor mendy & Richstond[T995i and references 
therein), but a much improved stellar-dynamical m e asure- 
me nt was recently obtain ed by IGebhardt & Thomas! (120091) 
and IGebhardt et all (1201 ll) using three-integral, axisymmetric 
Schwarzschild models that also include t he ma ss distribution 
of the dark halo. More specifically, IGebhardt et al.l (120111) 
applied orbit-based models to high angular resolution spec- 
troscopy obtained with Gemini's Near-Infrared Integral Field 
Spectrome ter assisted by a daptive o ptics and addition al large- 
scale data dEmsellem et al . 2004; M urphvet al.l2011l) . finding 
Mbh = (6.6 ± 0.4) X WMq. 

In addition to the stellar-dynamical Mbh measurements, 
there are also two gas-dynamical measurements for the black 
hole in M87. In fact, M87 was the first gas-dynamical tar- 
get with the Hubble Space Telescope (HST). Using spec- 
tra acquired with the Faint Object Spectrograph (FOS), 
[Harms et al.] (119941) modeled the emission-line kinematics as 
a rotating disk to i nfer a central mass of (2 .4 ± 0.7) x 10^ M©. 
A few years later, i Macchetto et al.l (11997 1) used spectroscopy 
from the Faint Object Camera (FOC) and rotating, thin-disk 
models that, for the first time, incorporated the important ef- 
fects of the propagation of light through the telescope and 
spectrograph opti cs to determin e Mb h = (3.2 ±0.9) x 10^ M©. 

The work of iHarms et all ([T994) and IMacchetto et al.1 
( 119971) were important milestones for HST and the field of su- 
permassive black hole detection. However, neither study was 
able to fully map out the kinematic structure of the gas disk. 
[Harms et al. ( 1994) relied on just six spectra obtained with the 
FOS and Macchetto et al. (1997) used FOC long-slit obser- 
vations at three, non-contiguous, parallel positions that only 
covered a fraction of the disk. Also, the disk inclination has 
been a sou rce of uncertainty in the M87 gas-dynamical black 
hole mass. "Harms et al.'('1994') adopted a value of / = 42° ±5°, 
which was measured by Ford et al. ( 1994) from HST Ha nar- 
rowband images, while IMacchetto et al.l (1 1997!) allowed the 



inclination to be a free parameter in their dynamical models, 
but found a wide range of acceptable values with 47° < / < 
65°. Perhaps more importantly, later investigations of cir- 
cumnuclear gas disks revea led the prevalence of an intrin- 
sic ve locity dispersion (e.g., van der Marel & van den Bosch] 
119981; Verdoes Kleijn et al. 2000; Dalla Bonta et al. 2009). 
Although the origin of the intrinsic velocity dispersion is not 
known, it may be dynamically important, contributing pres- 
sure support to the disk. Models that treat such systems as 
simple, thin, dynamically cold disks will underestimate Mbh 
(Barth et al. 2001; Neumayeret al. 2007; Walsh et al. 2010). 
Only by carrying out detailed modeling of the emission-line 
profiles, including the contributions to the line widths from 
rotational and instrumental broadening, will it be possible to 
determine if there is significant intrinsic velocity dispersion in 
the M87 disk and its subsequent effect on the gas-dynamical 
Mbh measurement. 

Given that M87 is located at the extreme high-mass end of 
the Mbh - host galaxy relations and is a very important test 
case for a gas/stars comparison study, the factor of ~2 dis- 
crepancy between the best current stellar and gas-dynamical 
masses is troubling. Furthermore, the mass of the black hole 
in M87 is a crucial parameter for making inferences about 
the size of the innermost stable circular orbit of the black 
hole and the black s pin from Very Long B aseline Interfer- 
ometry observations (iDoeleman et al.l [20121) . In this paper, 
we revisit the gas-dynamical black hole mass measurement. 
We use new HST data acquired with the Space Telescope 
Imaging Spectrograph (STIS) and carry out more compre- 
hensive gas-dynamical modeling than has previously been at- 
tempted for this galaxy. We describe the observations in ^ 
the emission-line kinematics in ^ the modeling procedures 
in ^and[5j the results and error analysis in §[6] and conclude 
in ^ Throughout the paper, we assume a distance to M87 
of D = 17.9 Mpc in order to allow for a dir ect comparison 
to the stellar-dynamical Mbh measurement bv lGebhardt et al.l 
([20Tlh . 

2. HST OBSERVATIONS AND DATA REDUCTION 

We obtained new HST STIS observations of M87 on 
2011 June 7. The observations have a wavelength scale 
of 0.554 A pixel"' and a spatial scale of 0"0507 pixel"'. 
We used the G750M grating, centered on 6581 A, to pro- 
vide coverage of the Ha spectral region, which includes 
the [O I] AA6300,6364, Ha, [N II] AA6548,6583, and [S II] 
AA6716,6731 emission lines. However, the [O I] and [S II] 
lines were weak, except at locations very close to the nucleus, 
and we focus solely on the Ha and [N II] lines for the re- 
mainder of the paper We used the El aperture position to 
place the nucleus near the readout end of the CCD, to avoid 
charge transfer efficiency losses. The STIS 52x0.1 slit was 
placed at five parallel positions without any space between 
adjacent positions, and was oriented at a position angle (PA) 
of 51° [with in 15° of the nuclear ga s disk's major axis as de- 
termined by [M^chettoeral]([T997[)]- Between 3-5 dithered 
exposures were acquired at each slit position to aid in the re- 
moval of cosmic rays, leading to exposure times that ranged 
from 1521 -291 1 s per slit position. 

The data from each subexposure were reduced by trim- 
ming the overscan region, subtracting the bias and dark files, 
and applying flat-field correc tions. After this init ial process- 
ing, we used LA- COSMIC dvan DokkumI 1200 Ih to remove 
cosmic rays and hot pixels. We completed the reduction of 
each subexposure by wavelength and flux calibrating the data 




Fig. 1.— HST WFPC2/PC F547M (left) and continuum-subtracted 
Hq-i-[N II] (right) images of the M87 nucleus. The images have been ro- 
tated such that the STIS instrumental y-axis points upward, and each box is 
1''7 on a side. The five STIS slit positions are overplotted on the continuum- 
subtracted image, and the length of the rectangles show the region over which 
the emission lines could be measured. 



and rectifying for geometric distortions. The reductions were 
carried out using individual IRAFO tasks within the standard 
Space Telescope Science Institute (STScI) pipeline, with the 
exception of the cosmic-ray cleaning step. The reduced, geo- 
metrically rectified, subexposures were aligned and combined 
using the IRAF tasks IMSHIFT and IMCOMBINE in order to 
produce the final two-dimensional (2D) STIS image at each 
of the five slit positions. 

Furthermore, from the HST archive we retrieved Wide 
Field Planetary Camera 2 (WFPC2) F658N, F547M, and 
F702W images, which were originally made under programs 
GO-5122 and GO-5476 and had the M87 nucleus placed 
on the PC detector The images in each filter were taken 
as a sequence of two individual exposures with total ex- 
posure times of 2700s, 800 s, and 280 s for the F658N, 
F547M, and F702W images, respectively. The subexposures 
were combined and the cosmic rays were rejected using the 
IMCOMBINE and LA-COSMIC IRAF tasks. We then cre- 
ated a continuum-subtracted Ha-i-[N II] image by rotating the 
F658N and F702W images to a common orientation, and sub- 
tracting a scaled version of the F702W image from the F658N 
image. The scale factor was chosen such that a background 
region in the continuum-subtracted image would have a mean 
flux as close to zero as possible. In Figure [T] we present the 
emission-line image with the location of the STIS slits over- 
laid, as well as the F547M WFPC2 image. A compact central 
source and an extended disk structure with several spiral-like 
wisps is seen in t he contin uum-subtracted image, as was pre- 
viously found bv lFord et al.l (II9941) . 

3. EMISSION-LINE KINEMATICS 

From each of the 2D STIS images, we extracted spectra 
from individual rows out to ~0"5 from the slit center, which 
is as far as the emission lines were detectable. At locations 
where the emission lines became weak, we binned together 
three STIS rows in order to improve the signal-to-noise ratio 
(S/N). Prior to measuring the emission-line kinematics, we 
subtracted the starlight continuum on a row-by-row basis by 
fitting a line between wavelengths 6295 - 6862 A excluding 
the regions near the emission lines, then subtracting this fit 
from the spectrum. Due to the low S/N in the continuum and 
the small wavelength range, a straight line is a sufficient de- 
scription of the continuum near Ha and we do not attempt a 

* IRAF is distributed by the National Optical Astronomy Observatory, 
which is operated by the Association of Universities for Research in Astron- 
omy under cooperative agreement with the National Science Foundation 



more complicated subtraction. Below we detail our spectral 
fitting procedure and describe the observed velocity fields. 

3.1. Spectral Fitting Method 

We simultaneously fit three single Gaussians to the 
continuum-subtracted Ha and [N II] emission lines using the 
MPFIT library in IDL (Mai-kwardt 2009), which employs a 
Levenberg-Marquardt least-squares minimization technique. 
We then performed a Monte Carlo, adding random Gaussian 
noise to the spectrum based upon the model residuals, and 
refit the new spectrum to determine the velocity, velocity dis- 
persion, and flux for each of the three emission lines during 
each iteration. From the resulting distributions, we were able 
to determine whether the [N II] A6548 and A6583 velocities 
were consistent within 2a. If the velocities were consistent, 
we calculated the mean and the 2a errors from the combined 
distribution of the [N II] velocities, which was taken to be the 
final velocity and error for that spectral row. The same proce ss 
was adopted for the [N II] velocity dispersions and for 3 x the 
flux of the [N II] A6548 line and the flux of the [N II] A6583 
line. If instead the kinematics of the [N II] lines were not con- 
sistent, we fit progressively more constrained models (such as 
tying together the [N II] velocity widths) until the [N II] kine- 
matics were consistent. With our approach, we are assuming 
that the [N II] velocities should be equal, the [N II] velocity 
dispersions should be the same, the fluxes of the [N II] lines 
should maintain a 3 : 1 ratio (as dictated by atomic physics), 
and that any differences that are seen are due to fitting error. 
We chose not to fit Gaussian models that enforce these con- 
straints from the start because we found doing so produced 
very small errors on the kinematics - typically we measured 
velocity errors of only a few km s"' for the good S/N spectra 
(while the STIS pixel scale is ^25 km s"'). 

We also tested whether the spectral fits could be signifi- 
cantly improved with the addition of a broad Ha component 
described by a single Gaussian. Including a broad component 
either caused the x^ per degree of freedom (xj,) to decrease 
by less than 15%, produced a broad component with zero flux, 
or resulted in unreasonably small fluxes and/or line widths 
for the Ha and [N II] narrow components. We cannot fully 
rule out the existence of a compact broad-line region, but the 
data do not require the presence of a physically distinct broad- 
line component, s imilar to the giant elliptical galaxy M84 
(I Walsh et al.ll20IOl) . In both these galaxies, the very broad 
central line widths are consistent with coming from the inner 
portion of the rotating disk. Thus, the kinematics presented 
here come from single Gaussian fits to the Ha and [N II] lines 
only without the inclusion of a model for a broad component. 

While the above fitting method worked well for many of 
the observed spectra, the spectra extracted near the center of 
slit positions 2, 3, and 4 were complex with severely blended 
emission lines. In addition, some of these rows showed 
artifacts, in the form of a periodic fluctuation in the flux, 
from the geometric rectification step of the data reduction 
jDressel et al.ll2007h . Moreover, spectra extracted near the 
center of slit position 5 exhibited blended [N II] A6548 and 
Ha emission lines and a noisy [N II] A6583 line that may 
be double-peaked. It was therefore difficult to decompose 
the spectra into Ha and [N II] components with accurate and 
unique mean velocities, velocity dispersions, and fluxes, even 
when applying constraints to the three single Gaussians that 
were being fit. Consequently, we consider the measurements 
from the central rows of slits 2, 3, 4, and 5 to be unreliable and 
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Fig. 2. — Example fits to the Hq and [N II] emission lines. Blue dotted 
lines show the three single Gaussians that were fit, and the red solid line is 
the sum. The top two panels depict the range in quality of spectra from which 
kinematic measurements were made. Close to the nucleus and near the center 
of slit position 5. it was difficult to accurately and uniquely decompose the 
spectra into into the Ha and [N II] components, as can be seen in the bottom 
two panels. Therefore, the kinematics measured from these central rows were 
excluded during the dynamical modeling. 

we do not use them to constrain our gas-dynamical models. 
The innermost kinematic measurements provide relatively lit- 
tle useful information when constraining the disk models be- 
cause the line profiles near the nucleus are dominated by ro- 
tational broadening and blurring by the point-spread function 
(PSF). For completeness, however, we incorporate into our fi- 
nal error budget the effect on the black hole mass when the 
measurements deemed as unreliable are used to constrain the 
dynamical model as well. In Figure |2] we show example fits 
to spectra extracted from slit positions 3 and 5, including rows 
where the kinematic measurements are considered unreliable. 

3.2. Observed Velocity Fields 

From the spectral fitting, we were able to map out the 
emission-line velocity, velocity dispersion, and flux as a func- 
tion of location along the slit for each of the five slit posi- 
tions. The velocity curves reveal that the gas participates in 
regular rotation, and the gas to the northeast of the nucleus is 
redshifted relative to the galaxy's systemic velocity while the 
gas located to the southwest is blueshifted. There is a steep 



velocity gradient of ^1200 km s"' across the inner Q'.'l of 
slit position 3. The observed velocity dispersions, which in- 
clude contributions from rotational and instrumental broaden- 
ing, range from ~'150 - 450 km s"', although larger velocity 
dispersions are likely present as the line widths continue to 
rise toward the nucleus. However, we were unable to obtain 
reliable measurements of the central velocity dispersions be- 
cause the emission lines were severely blended together in the 
nuclear spectra. In addition, the observed velocity dispersion 
profile appears asymmetric, perhaps as a result of interactions 
with the one-sided jet in M87. For example, the gas has ve- 
locity dispersions of 461 km s"' and 263 km s"' at locations 
-0."1 and 0."1 away from the nucleus in slit position 3. 

4. GAS-DYNAMICAL MODELING: OVERVIEW 

The observed velocity fields are modeled assuming that the 
gas participates in circular rotation in a thin disk. We provide 
a brief overview of the modeling here, but refer the reader to 
iBarfli et al.l (1200 ll) and I Walsh et al.l (12010 ) for further details. 
At each radius (r) in the disk, the circular velocity is calcu- 
lated relative to the galaxy's systemic velocity (usys) based on 
the enclosed mass, which depends on the black hole mass, 
the stellar mass profile, and the mass-to-light ratio (T). The 
model velocity field is created on a pixel grid that is oversam- 
pled relative to the size of a STIS pixel, and is projected onto 
the plane of the sky for a given value of the disk inclination 
angle (/). 

The intrinsic line-of-sight velocity profiles are assumed to 
be Gaussian before passing through the telescope optics and 
are generated on a velocity grid with a bin size that matches 
the STIS scale. The intrinsic line profiles are centered on the 
projected line-of-sight velocity at each point on the model grid 
and weighted by the emission-line flux distribution. In addi- 
tion, the widths of the intrinsic line profiles are set based upon 
contributions from the thermal velocity dispersion of the gas 
and the velocity dispersion due to turbulent motion internal to 
the gas disk. 

Next, the model velocity field is "observed" in a man- 
ner that matches the STIS observations. This synthetic ob- 
servation includes accounting for the blurring by the tele- 
scope PSF by convolving a Tiny Tim ( Krist & Hook|[2004h 
PSF model with each velocity slice of the line profile grid. 
Also, we propagate the line profiles through the finite width 
of the STIS slit by applying shifts to the velocity depending 
on the location along the slit width that the photon enters 
dMacieiewski & Binne"vll2001h . After rebinning to the STIS 
pixel size, we are left with a model 2D spectral image simi- 
lar to the observed STIS data. At this stage, we convolve the 
model 2D STIS image with the Tiny Tim CCD charge diffu- 
sion kernel in order to account for the charge that is spread 
between non-subsampled, neighboring pixels. 

Finally, we extract spectra on a row-by-row basis from the 
model 2D image and fit a single Gaussian to the emission line, 
analogous to the measurements that were made from the data. 
We thus are able to measure the model velocity, velocity dis- 
persion, and flux as a function of position along the slit. The 
best model is taken to be the one that most closely matches 
the observations in the x^ sense. We use the downhill simplex 
algorithm by Press et all (119921) to minimize x^ and determine 
the best-fit parameters. We separately optimize the fits to the 
observed line widths and fluxes before determining the black 
hole mass from the model fit to the observed radial velocities. 
When fitting to the observed velocities, the model parameters 
are Mbh, T, /, Ws^s, d, ^Coffset, and ^offset- Here, 9 is the angle 



between the slit and the projected major axis of the disk, and 
•^offset and joffset are the offsets between the center slit and the 
nucleus. 

5. GAS-DYNAMICAL MODELING: APPLICATION TO M87 

Gas-dynamical models for M87 were calculated on a pixel 
grid oversampled by a factor of ,s = 6 relative to the STIS pixel 
size. We found that preliminary models with a subsampling 
factor of .? = 6 were able to best reproduce the observed veloci- 
ties and emission-line fluxes, and could be run in a reasonable 
amount of time. Thus, we selected s = 6 for the final models. 
We also incorporated the effect on Mbh of our subsampling 
fac tor c hoice into the error budget, which will be summarized 
in ^6.11 In addition, we used a PSF generated with Tiny Tim 
for a monochromic filter passband at 6600 A. The model is an 
s = 6, 0"4x0."4 portion of the fufl STIS PSF. Such a model 
excludes part of the extended PSF wings, with 12% of the to- 
tal flux being omitted. As will be shown in 36.11 the black 
hole mass is not strongly dependent on the size of the PSF 
model, but because the computation time increases rapidly 
with increasing PSF size, we chose to use a 0."4-diameter 
PSF. Finally, we fit disk models over the entire radial range 
covered by the observed kinematics. With a fitting radius of 
0."5, not only is there an adequate number of measurements 
to constrain the dynamical models, but also the disk model 
continues to provide a good description of the observations 
out to this radius. As discussed previously, we exclude the 
uncertain kinematics measured from the central rows during 
the modeling. 

Below we provide additional details about the M87 gas- 
dynamical models, including the stellar contribution to the 
circular velocities, the parameterization of the intrinsic 
emission-line flux distribution, the components that make up 
the width of the line profiles, and our approach for allowing 
the turbulent velocity dispersion to be dynamically important. 

5.1. Stellar Mass Model 

We determined the stellar contribution to the gravitational 
poten tial using the V-b and surface brightness profile of M87 
from Kormendv et al.l (12009). The profile combines mea- 
surements from ground-based images with high-resolution 
HST WFPCl/PC data, which were originally presented in 
iLauer et aTl (119921) and have already been deconvolved by 
the telescope PSF. In addit ion to the surface brightness pro- 
file, iKormendy et alJ (1200 9) present the ellipticity (e) of the 
isophotes as a function of major axis radius. At radii within 
180", the isophotes are nearly circular, with e < 0.2. There- 
fore, we deprojected the surface brightness profile assum- 
ing that the intrinsic density distribution of the central re- 
gions could be modeled as the sum of spherically sym- 
metric components with Gaussian profiles. A similar ap- 
proach has been used to derive the luminous density distri- 
bution for past g as-dynamical black hole mass measurements 
(e.g.,[S arzi et aL lgOOlt lBar5rerani200TI: iCoccato et al.ll2006l: 
IWalsh e t al. lOW). To carry out the deprojection, we used 
the Mu lti-Gaussian Expansion (MGE) software developed by 
iCappellari (2002a), which is based on the methods described 
in Em sellem et al. (1994 ). In addition, the M87 nucleus is 
classifie d as a Type 2 L ow Ionization Nuclear Emission-line 
Region dHo etal.ll 19971) . We attribute the innermost compact 
Gaussian component of the MGE model, with a dispersion of 
3.5 pc, to light from the AGN, and exclude this component 
from the stellar mass distribution. 




Fig. 3 . — The cii'cular velocity as a function of radius due to tlie presence of 
a 3.0 X lO' Mq black hole (dot-dashed black curve) and the stars with T = 4 
(solid black Une) and T = 1 (grey solid line). The solid red curve shows the 
sum of the circular velocities from the black hole and the stars when assuming 
T = 4 (V-band, solar). The dashed vertical grey line indicates the largest 
radius probed by the STIS data. Over this range, the black hole dominates 
the gravitational potential and the stellar contribution is neghgible. 



The final product of the multi-Gaussian deprojection of the 
surface brightness profile is the stellar contribution to the cir- 
cular velocity in the galaxy, which is shown in Figure |3] The 
figure also provides the circular velocities due to the presence 
of a 3.0 X 10'' Mq black hole, which is a plausible lower limit 
on the black hole mass from our analysis. Clearly, the black 
hole's sphere of influence is well resolved by our STIS ob- 
servations, which is further supported by the nearly Keple- 
rian shape of the multi-slit velocity curves. In other words, 
the black hole dominates the galaxy's gravitational potential 
and the stellar contribution is negligible over the radial range 
probed by the STIS kinematics. Consequently, we found that 
our dynamical models are insensitive to the stellar mass-to- 
light ratio, and we chose to fix T = 4 (V-band, solar units) 
during the model fitting. This value was selected using the 
B-V color for the galax y from the Third Referenc e Cata- 
logue of Bright Galaxies ( de Vaucouleurs et al.lll991h and the 
relations presented in Bell et al. (2003). Such a stellar mass- 
to-light r atio is typical ly found at the centers of giant elliptical 
galaxies (ICappellari e t al. 2006). We also tested the effect on 
Mbh when another fixed value for T is adopted, as well as 
when T is allowed to be a free parameter. The results will be 
described in 36.11 

5.2. Emission-Line Flux 

The line-of-sight velocity profiles are weighted by the 
emission-line flux distribution at each point in the disk. Com- 
monly, the flux distribution is described by analytic func- 
tions, or by creating a continuum-subtracted Ha + [N II] im- 
age that has been deconvolved by the telescope PSF (e.g. , 
'Barth et al."2001t iPastorini et al.ll2007t iMarconi et al.l 120061: 
Wold et al. 2006!). We experimented with directly folding a 
deconvolved continuum-subtracted image constructed from 
WFPC2/PC F702W and F658N images into the modeling, as 
well as using ten analytic forms composed of two and three 
Gaussian [Fir) = Foexp(-r^/2rj)] and exponential [Fir) = 
Foexp(-r/rj)] components. We used analytic models that rep- 



TABLE 1 
Emission-Line Flux Parameters 



Component 


Fo 


(pc) 




Veen 

(") 


PA 

(°) 


b/a 


Gaussian 
Gaussian 


26.5 
LO 


6.6 

23.7 


0.05 
0.00 


-0.02 
0.02 


14 
191 


0.13 
0.77 



Note. — The amplitude, Fq, is in arbitrary units. The center of the 
ellipse is described by Xccn and Vccn relative to the location of the black 
hole. The position angle, PA, is in units of degrees clockwise, with PA = 
pointing along the length of the STIS slit. 

resented intrinsically circularly symmetric disks, which pro- 
duced concentric elliptical isophotes with constant position 
angle and axis ratio, and more complicated models in which 
the isophotes of the individual components were allowed to 
have different centers, position angles, and axis ratios. The 
parameters associated with each of the analytic functions were 
determined by calculating disk models following the methods 
described in ^ then optimizing the fit to minimize the x^ 
computed from the comparison of the model and observed 
emission-line fluxes. 

We found that using a deconvolved continuum-subtracted 
image and analytic models that represented intrinsically cir- 
cularly symmetric disks returned unacceptable fits to the 
observed emission-line flux and velocities, and we do not 
consider these descriptions further While deconvolved 
Ha + [N II] images have been folded i nto gas-dynamical 
models successfully in the past (e.g., iBarth et aTl 12001 1: 
iDalla Bonta et alj 120091) . we experienced difficulties during 
the deconvolution process, which was performed using the 
IRAF task LUCY. The deconvolution produced an artifact in 
the form of a hole near the nucleus and subsequently the 
model fluxes were smaller than the observed fluxes for a por- 
tion of the disk. The remaining five analytic models, how- 
ever, were able to match the observed flux distribution, and 
we searched for the simplest model that would simultaneously 
produce model velocities closest to the observed radial veloci- 
ties. Ultimately, we determined that the best parameterization 
of the intrinsic emission-line flux was the sum of two Gaus- 
sian components. In Tablefl] we present the amplitudes (Fo), 
scale radii (rj along the major axis, position angles, axis ra- 
tios (b/a), and centroid positions for each of the two compo- 
nents. We also include the effect of using the other four an- 
alytic flu x mo dels in the final error budget for the black hole 
mass in ^6.11 

5.3. Line Widths 

We initially ran models in which the widths of the emission- 
line profiles originate from rotational, instrumental, and ther- 
mal broadening. The rotational and instrumental broadening 
effects are explicitly included in our dynamical models when 
the line profiles are propagated through the telescope and in- 
strument. In particular, rotational broadening occurs because 
light from different portions of the disk get blended together 
within the slit, while instrumental broadening occurs because 
of PSF effects, the charge diffusion between neighboring pix- 
els, and the velocity shifts that result from the finite slit width. 
Prior to propagating the line profiles through the telescope op- 
tics, we also include a thermal contribution to the velocity dis- 
persion of (Tth = 10 km s~', which is expected for gas with a 
temperature of salO'* K. 

In some parts of the disk, the predictions from these 
initial models were smaller than the observed veloc- 



ity dispersions. Similar results have been found when 
constructin g gas-dynamical models of other galaxies as 
well (e .g., iv an der Marel & van den Bosch'^ggS"; 'Barth et alJ 
200U rV^ oes Kleiin et al. 2002; Tadhunter et al. 2001 



Dalla Bonta et al . 2009), although in some instances the ro- 



tational and instrumental broadening alone is enough to ex- 
plain the observed line widths (e.g., [M arconi et al. 120011; 
Capetti et al.l l2005t [Atkinson et all 120058: ,de Francesco et alj 
2006, 2008). Therefore, we added a projected intrinsic ve- 
locity dispersion (dp) in quadrature to (Tth before propagat- 
ing the line profiles through the telescope and spectrograph. 
The projected intrinsic velocity dispersion was determined by 
constructing disk models as discussed in ^and fitting to the 
observed [N II] line widths. With our early model runs, we 
tried characterizing Up as a constant + exponential of the form 
(To + (Ti exp(-r/r()), which has often worked well for other 
galaxies, but found that best-fit value for a\ was ^0 km s"'. 
Thus, in the final model, the intrinsic velocity dispersion was 
parameterized as a constant over the disk, with ijp = 170 km 



5.4. Asymmetric Drift Correction 

Currently, there is no clear consensus for how to inter- 
pret the physical nature of the intrinsic velocity dispersion, 
hence we explore two scenarios that should bracket the range 
of possible masses for the black hole in M87. In one case, 
we assume that the intrinsic velocity dispersion does not 
affect the circular velocity, which could result if the in- 
trinsic velocity dispersion is due to local microturbulence 
but the bulk motion of th e gas remains in circular motion 
dvan der Marel & van den Bosch 1998). In the second case, 
we assume that the internal velocity dispersion is due to turbu- 
lent motion that supplies pressure support to the disk, causing 
the observed mean rotational speed (v^f,) to be smaller than 
the l ocal circular velocity (Vr) for a given black hole mass 
(e^ lCoccato et alj|2006t iNeumaver et al.ll2007l: IWalsh et al.1 
120101) . For M87, we estimate an upper bound to the black 
hol e mass by applyin g an asymmetric drift correction follow- 
ing lBarth et al.l(l2001h . Under the assumption that the gas mo- 
tions are isotropic in the r and z (cylindrical) coordinates, the 
correction is 



d\nv d\\\(j,. 



d\nr d\nr 



(1) 



We take the intrinsic emission-line flux distribution to be a 
proxy for the number density of tracer particles, or clouds in 
the gas disk, v. Since the asymmetric drift correction depends 
on an azimuthally-averaged value for the radial derivative of 
v, we refit the observed fluxes with an analytic function com- 
posedof two Gaussians (similar to the flux model presented 
in ^5.2b . but required that the isophotes have the same cen- 
ter, position angle, and axis ratio. For a given intrinsic radial 
velocity dispersion, ar, we compu ted a^, and the p rojected 
velocity dispersion as outlined by Barth et aP (120011) . before 
adding the projected velocity dispersion in quadrature to crth 
and propagating the line profiles through the telescope and in- 
strument. Although we initially attempted to parameterize a,- 
as a constant + exponential function, we found that the intrin- 
sic radial velocity dispersion was best described as constant. 
By calculating disk models with an asymmetric drift correc- 
tion, and fitting to the observed line widths, we found a best-fit 
value of cr,. = 191 km s"' . 



TABLE 2 
Disk Model Parameters 



(Mq) 



n 



Vsys 

(kms"') 



(°) 



^offset 



(") 



O.S!*;^) X lO' 42+5 



1335;3 



6±2 



-0.02 ±0.01 0.03:»<j;]2 



Note. — These results give the best-fitting parameters for models cal- 
culated without an asymmetric drift correction. The errors given for the 
black hole mass are Icr uncertainties, and the errors provided for the other 
parameters correspond to the range of values found from all the preUmi- 
nary models that were run. During the modeling, T was frozen at 4 {V- 
band solar units). A relative angle of 9 = 6° corresponds to a disk major 
axis position angle 45° east of north. 

6. MODELING RESULTS 

After optimizing the model fits to the observed emission- 
line fluxes and velocity dispersions, we fit disk models to the 
observed radial velocities. The results returned by the x^- 
minimization routine seemed to strongly depend on the ini- 
tial guesses for Mbh and / in particular. Therefore, to ensure 
a global minimum was found, we used a two-step process. 
First, we ran a coarse grid of models with fixed Mbh and / 
values while the other parameters (9, Ujys, ^offset, ^offset) were 
allowed to vary. Second, we used the parameter values of the 
top five models (with the lowest x^) from the grid search as 
the starting points to do five runs in which all the parameters 
were allowed to vary. By doing these five runs, we were able 
to distinguish between local and global minima. 

The results of the best-fit model without an asymmetric drift 
correction are summarized in Tableland the model predic- 
tions are shown alongside the data in Figures |4] and |5] We 
measure a black hole mass of Mbh = 3.5 x 10^ Mq. The 
other model parameters, including the disk inclination an- 
gle, are well constrained. During the numerous preliminary 
model runs, the parameters varied at most by 7°, 9 km s~', 
2°, 0"01, and 0"04 for /, Usys, 9, JCoffset, and ^offset, respec- 
tively. The model has x^ = 66, found by comparing the model 
predictions to the observed velocities. Given that 62 veloc- 
ity measurements were used to constrain 6 model parameters, 
the x^ P6r degree of freedom is xl = 1-2- As can be seen 
in Figure H] the model is able to reproduce the shape of the 
observed velocity curves and the emission-line flux profiles 
very well. At some locations in the disk, an intrinsic velocity 
dispersion is needed to match the data, although this velocity 
dispersion is relatively small given the large velocities found 
in M87. Moreover, our model for the intrinsic velocity dis- 
persion still severely underestimates the velocity dispersions 
located -0."10 and -0"15 away from the center of slits 3 and 
4. We will examine this discrepan cy an d the possible impact 
on the inferred black hole mass in ^6.11 

6. 1. Uncertainty in the Black Hole Mass 

We measured the statistical errors on Mbh by holding Mbh 
fixed while allowing the other parameters to vary such that x^ 
was minimized. During this process, we froze Mbh to val- 
ues between (1.0-9.7) x lO'' Mq in steps of 0.1 x 10^ Mq. 
We searched for the range of Mbh values that caused the min- 
imum x^ to increase by 1.0 and 9.0 to estimate the la and 
3(7 model fitting uncertainties. Figure |6] presents the results 
of these disk models, specifically the x^ variation with Mbh 
after marginalizing over the other parameters. Ultimately, we 
determine that the range of black hole masses is (3.4-4.0) x 
lO"* Mq (la uncertainties) and (3.0-6.8) x lO** Mq (3cr un- 



certainties). 

These Mbh errors account for the random noise in the gas 
dynamical models and the uncertainties associated with the 
parameters /, Usys, 9, Xoffset, and ^offset- However, there are ad- 
ditional sources of uncertainty that are not incorporated into 
the statistical errors. Below, we explore these sources of un- 
certainty and their impact on Mbh, in order to arrive at the 
final range of possible black hole masses forM87. 

Subsampling Factor. We ran gas-dynamical models on pixel 
grids that were subsampled by a factor of i = 2- 10, in in- 
tervals of s = 2. We found that the black hole mass varied 
between (3.5 - 3.9) x 10^ Mq with a root-mean-square (rms) 
scatter of2. 3 x 10** M©, which is 7% of the best-fit black hole 
mass. 

PSF Size: When calculating dynamical models using a Tiny 
Tim PSF with diameters of 0."3, 0."4, 0."5, and 0"6, the rms 
scatter in the black hole mass was 1.2 x 10** Mq. The PSF 
size has a minimal impact on the inferred black hole mass, as 
the scatter is just 3% of the best-fit mass. 

Stellar Mass-to-Light Ratio: The stellar contribution to the 
galaxy's gravitational potential is negligible over the region 
probed by the STIS data ('^40 pc). As a result, our models 
were unable to constrain the stellar mass-to-light ratio, and 
preliminary model runs returned a wide range of values for 
T. We chose to fix T = 4 (V-band, solar) during the final 
model run, but also tested freezing T to 9, which is the stel- 
lar mass-to-light ratio found for M87 b y Murphy et al. (201 Ij) 
and adopted by Gebhardt et al.l (1201 Ih . and allowing the pa- 
rameter to float. The test produced in an rms scatter in the 
black hole mass of 1 .3 x 10^ Mq, or 4% of the best-fit Mbh- 

Emision-line Flux Model: The intrinsic emission-line flux 
distribution was modeled as the sum of two Gaussian com- 
ponents, whose isophotes had different centers, position an- 
gles, and axis ratios. An additional four analytic forms of 
varying complexity were also able to reproduce the observed 
emission-line fluxes: two exponentials, one Gaussian + one 
exponential, three Gaussians, three exponentials, and two 
Gaussians + one exponential. These analytic functions also 
represented elliptical isophotes with different centers, posi- 
tion angles, and axis ratios. Although the quality of the fit 
to the observed velocities varied significantly depending on 
the intrinsic flux distribution used (x^ ranged from 66 - 130), 
there was a small effect on the black hole mass. The disk mod- 
els with different parameterizations of the emission-line flux 
distribution returned best-fit masses between 3.1 x lO' Mq 
and 4.2 x lO' Mq, with an an rms scatter of 5 .0 x 10^ Mq, or 
7% of the best-fit black hole mass. 

Differences Between the [N II] and Ha Kinematics: Gas- 
dynamical models were run using the [N II] kinematics be- 
cause [N II] A6583 was the strongest emission line in our spec- 
tra. If gas-dynamical models are instead fit to the Ha kinemat- 
ics, the black hole mass increases by 1 1% to 3.9 x lO' Mq. 

Including all Kinematic Measurements: We had difficulty fit- 
ting the spectra extracted near the center of slits 2-5. We con- 
sider the kinematics and emission-line fluxes measured from 
these rows unreliable and did not use them to constrain our 
final disk model. When these velocities, velocity dispersions, 
and emission-line fluxes are instead included in the fit, the 
best-fit black hole mass decreased by 1 1 % to 3 . 1 x 10^ Mq . 

Asymmetric Drift Correction: By assigning a dynamical ori- 
gin to the intrinsic velocity dispersion and calculating an 
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Fig. 4. — The emission-line velocity, velocity dispersion, and flux (black points) measured as a function of location along the slit for the five slit positions. The 
numbers at the top left of each plot correspond to the slit locations shown in Figure [T] while the negative and positive Y Offsets con'espond to the bottom and 
top half of the slits, respectively, depicted in Figure [T] We defined Y Offset = to be the spectral row with the largest continuum flux. Grey points indicate the 
velocity measurements that were considered unreliable and are plotted for reference, while the corresponding velocity dispersion and flux measurements are left 
off the figure. Spectral rows that were binned together to improve the S/N are shown with en'or bars along the A-axis. The spectra from the bottom portion of 
slit 1 had very poor S/N and we were unable to obtain kinematic measurements. The best-fit dynamical model for M87 without an asymmetric drift correction 
is overplotted with the red solid line, and the best-fit systematic velocity is displayed with the dashed black fine. A single scahng factor has been applied to the 
model emission-Hne fluxes such that the median values of the model and data would match. The dotted red line shows the model line widths without the inclusion 
of an intrinsic velocity dispersion. 




Fig. 5. — Continuum-subtracted Ha 4- [N II] region of the STIS data (left) and the synthetic 2D spectrum of the [N II] A6583 emission line from the best-fit 
dynamical model without an asymmetric drift correction (right) at each of the five slit positions. Each box is 93 A along the dispersion (horizontal) direction and 
2f'2 along the spatial (vertical) direction. The numbers correspond to the slit positions in Figure[T] 
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Fig. 6. — Results of the disk models run without an asymmetric drift cor- 
rection by fixing Mbh while allowing the other paraineters to vary such that 
X^ was miniinized. The dashed and dotted lines denote where the minimum 
X^ has increased by 1.0 and 9.0, respectively, which corresponds to the Icr 
and 3(T uncertainties for one degree of freedom. 



asymmetric drift correction as discussed in § 15.41 the black 
hole mass increased by 6% to 3.7 x 10^ Mq. Such a minor 
change to Mbh is not surprising considering the small amount 
of intrinsic velocity dispersion relative to the large velocities 
that are seen in M87. Compared to the dynamically cold, thin- 
disk model, the model with an asymmetric drift correction 
was a worse fit to the observed velocities, with a x^ = 71. 

The predictions from the model with asymmetric drift cor- 
rection, however, continue to be significantly smaller (by 
^140 - 240 km s"') than the observed dispersions located 
-0"10 and -0"15 away from the center of slits 3 and 4. In 
order to establish the maximum possible increase in the black 
hole mass, we searched for a new parameterization of the in- 
trinsic velocity dispersion that would match these few points, 
at the expense of being able to reproduce the observed veloc- 
ity dispersions at the other locations. Since the distribution of 
the observed line widths is asymmetric, and the velocity dis- 
persions at negative Y Offsets are generally larger than those 
measured at positive Y Offsets, we modified the input kine- 
matics by simply replacing the observed velocity dispersions 
from the positive Y Offset side with the values from the nega- 
tive Y Offset side. We then fit gas-dynamical models with an 
asymmetric drift correction first to the modified line widths 
and then to the velocities. We determined that the intrinsic 
velocity dispersion was best characterized as a constant + ex- 
ponential function, with cto = 202 km s"', ai = 2040 km s"', 
and ro = 4.4 pc, and that the black hole mass increased by 20% 
to Mbh = 4.2 x 10^ Mq. As can be seen in Figure|2l the new 
model is able to sufficiently match the high velocity disper- 
sions observed at Y Offsets of -0"10 and -0"15 from slits 3 
and 4. 

The formal model fitting error and the additional sources 
of uncertainty were added in quadrature to derive the final 
range of black hole masses. Thus, we find the M87 black hole 
mass is (2.8 - 4.4) xlO^ Mq (Ict uncertainties) and (2.7 - 
6.9) X 10*^ Mq (3cr uncertainties) with a best-fit value of 3.5 x 
10^ Mq. We propagated the 6% increase in the black hole 
mass resulting from the asymmetric drift correction into the 
upper uncertainty of Mbh only. Even when using an extreme 
characterization of the intrinsic velocity dispersion to better 
match the observed line widths at Y Offsets of -0."10 and 
-0"15 in slits 3 and 4, the black hole mass does not increase 



beyond the final Icr uncertainties. 

7. DISCUSSION AND CONCLUSIONS 

Using newly acquired STIS observations of M87 with the 
0"l-wide aperture at five adjacent parallel positions, we have 
measured the [N II] velocity, velocity dispersion, and flux 
well within the black hole's dynamical sphere of influence. 
Through gas-dynamical modeling of the STIS data, we have 
determined that the mass of the black hole is Mbh = (3 -Sil^o 7) x 
10^ Mq. This black hole mass is similar to the expecta- 
tions fromtheMBH - bulge relations. With tr* = 324 km 
s"' (iGebhardt et al. 2011) an d a V-band l uminosity calcu- 
lated from My = -22.71 ( Lauer et al.ll2007h . the most recent 
calibrations of the Mbh ~ o"* relation suggest that Mbh ^ 
3 x 10^ Mq dCiraham & ScotdlMl iMcConnell &"Mall20Tl 
and the MBH-^bui correlation predicts Mbh = 2 x 10^ Mq 
(McConnefl & Ma 20T3|). 

Unl ike t he past gas-dynamica l studies of iHarms et al.l 
(|I99 4 and lMacchetto et al.1 (Il997h . we have mapped out the 
complete kinematic structure of the emission-line disk. We 
have constrained the gas disk inclination angle, which has pre- 
viously been a source of uncertainty, to be between 35° < / < 
47° with a best-fit value of / = 42°. In addition, our models 
included a detailed treatment of the effects of the telescope 
and instrument. For the first time, we have found that a rel- 
atively small amount of additional velocity dispersion inter- 
nal to t he gas disk is needed to reproduced the observed line 
widths. iMacchetto et alj (1 1 9971) did not find evidence for such 
an intrinsic velocity dispersion in the M87 disk, but the instru- 
mental resolution of their FOC observations was low (FWHM 
~ 430 km s"' ). Thus, the order-of-magnitude improvement in 
spectral resolution of STIS compared to the older FOC data 
has given us the sensitivity to detect an intrinsic velocity dis- 
persion. We considered the possibility that the intrinsic ve- 
locity dispersion provides dynamical support to the gas disk, 
and determined that inferred mass increases by just 6%. We 
incorporated this effect, along with a number of other sources 
of uncertainty, into the error budget. 

Before comparing our mass measurement to the prior mass 
determinations, there are some additional systematics to keep 
in mind when interpreting the results. Our thin-disk model 
matches the overall shape of the observed velocity curves 
very well and returns xi = ^ -2, however the data is not ran- 
domly scattered about the best-fit model at all locations. In- 
stead the model systematically deviates from the observa- 
tions in certain regions, at Y Offset = 0"25 - 0"4 at slit 
3, for example. Consequently, there may be some regions 
that depart from circular rotation, whose small-scale veloc- 
ity structure cannot be matched with any disk model. This 
is a common problem for gas-dynamical models in general, 
and makes estimating the uncertainties in the black hole 
mass challenging. Past studies have opted to rescale the 
error bars on the kinematics, such that the best-fit model 
would have a x^ ~ I , prior to measuring the statistical er- 
rors on Mrh (e.g iBarth et al.' '200!1, 'Atkinson et all 120051; 
Ide Francesco et alJ 12006: Dalla Bonta et al. 2009). For this 
work, we do not adjust the kinematic errors because our best- 
fit model already has x?- ~ 1- Despite the reasonable xi of 
our best-fit model, we note that the model is not a perfect de- 
scription of the gas kinematics everywhere in the disk. 

Furthermore, we treat the (small amount of ) intrinsic veloc- 
ity as dynamically important by applying an asymmetric drift 
correction. The correction is applicable in the limit of col- 
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Fig . 7 . — The predictions of a disk model ran with an asymmetric drift correction are compared to the data from slit positions 3 and 4; see Figure|4]for additional 
description. Here, an extreme intrinsic velocity dispersion has been chosen so that the model would match the large observed line widths located —Of 10 and 
— 0''15 away from the center of slits 3 and 4, at the expense of being able to reproduce the smaller observed hne widths at some of the other locations. This 
exercise provides an indication of the maximum possible increase in Mbh when assigning a dynamical origin to the intrinsic velocity dispersion. 



lisionless particles and for a^/vc ^ 1 . Here, we applied the 
correction to gas clouds, which are not collisionless, and to a 
system that has cr,./vc < 0.3 over the radial extent of the STIS 
kinematics. Also, the asymmetric drift correction depends 
upon the number density of clouds in the disk, v{r), which is 
unknown. While we assumed that the intrinsic emission-line 
flux distribution can be used as a proxy for v{r), this is at best 
a rough approximation. With these considerations in mind, 
our model with an asymmetric drift correction only provides 
estimation, albeit a useful one, of the dynamical influence of 
the intrinsic velocity dispersion. 

Our mass measurement is consistent with the two previous 
gas-dynamical mea sure ments of the black hole i n M87. The 
iHarms et all (119941) and iMacchetto et alj (119971) masses are 
Mbh = (2.9±0.8) x \(f Mq andMsH = (3.8± 1.1) x IO'^Mq, 
after scaling to our adopted distance. Of the remaining disk 
model parameters, our best-fit inclination angle is the same 
as the one used by Harm s et al.l (1199 4) and is consistent with 
the lower bound of accept able values (/ = 47° - 65°) found 
by IMacchetto et alj ( 119971) . while our best-fit systemic ve- 
locity (wsys = 1335 km s~') is larger than the ones used by 
IHarms et all (119941) (Usys = 1309 km s"') and determined by 
IMacchetto et al.l (119971) (.v^^^ = 1290 km s~'). For comparison, 
the recession velocities measured from optical lines given by 
the NASA Extragalactic Database span a wide range of val- 
ues, from 800- 1747 km s"'. We established that the pro- 
jected major axis of the gas disk is 45° east of north, which is 
roughly perpendicu lar to the jet and similar to position angle 
of 38° measured by Macchetto et al.l (119971) . 

The most recent M87 stellar-dynamical black hole mea- 
surement of Mbh = (6.6 ± 0.4) x 10'' Mq (Gebhai' dt et"al] 
120111) is a factor of two larger than our mass measurement, 
and there is a 2a discrepancy betw een our gas-dynamical 
mass and the iGebhardt et al.l (12011!) stellar-dynamical mass 
when accounting for the uncertainties associated with each 
measurement. Recent stellar-dynamical work has shown the 
importance in some galaxies of including a dark halo, in- 
corporating orbital libraries that better sample the phase - 
space occupied by tangential orbits, using triaxial geome- 
tries, and acc ounting for a spatially varying stellar mass- 
to-light rati o ("Gebhardt & Thomas' '2009^, 'Shen & Gebhardt' 
12010: Ivan den Bosch & de Zeeuw 20 10; Schulze & Gebhardt. 



201 I s iMcConnell et alllloTl . While IGebhardt et all dloTTb 
address the first two potential sources of systematic error, 
they do not examine the later two sources. In particular, 
they do not fit for the galaxy's intrinsic shape and orienta- 
tion, and instead use axisymmetric orbit-based models that as- 
sume an edge-on inclination. Therefore, it is still possible for 
the stellar-dynamical black hole mass measurement to change 
again if more general triaxial models are applied, although 
the magnitude and direction of the effect on the inferred black 
hole mass is unknown a priori. Also, Gebhardt et al. (2011]) 
construct stellar-dynamical models under the assumption of 
a constant stellar mass -to-light ratio, however, as shown by 
IMcConnell et alj (120 13h . spatial gradients in the mass-to-light 
ratio can have a non-negligible effect on Mbh- In contrast, 
our dynamical models rely on gas in circular orbits at small 
radii and are insensitive to the overall normalization of the 
stellar mass-to-light ratio ( ^6. IK making our gas-dynamical 
measurement immune to this uncertainty. 

Our gas-dyna mical measurement and the stellar-dynamical 
determination of IGebhardt et al.l (12011!) utilize high-quality, 
high-resolution observations and apply some of the most up- 
to-date modeling techniques, yet there continues to be a dis- 
crepancy in the M87 black hole mass. This highlights the 
need for carrying out more such cross-checks between the two 
methods. Currently, M87 is one of six consistency tests that 
have yielded a meaningful comparison between the gas and 
stellar-dynamical techniques, and in three cases (this study in- 
cluded) the stellar-dynamical mass exceeds the gas-dynamical 
determination by a factor of ^2 - 5. Since the stellar and 
gas-dynamical techniques are completely independent meth- 
ods with different systematic effects, only by applying both 
techniques to the same object can conclusions be made re- 
garding the consistency of the methods, the subsequent effects 
on the Mbh scaling relations, and the magnitude and distribu- 
tion of the intrinsic scatter of the Mbh - host galaxy relations. 
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